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Abstract 



For a freely evolving granular fluid, the buildup of spatial correlations in 
density and flow field is described using fluctuating hydrodynamics. The the- 
ory for incompressible flows is extended to the general, compressible case, 
including longitudinal velocity and density fluctuations, and yields qualita- 
tively different results for long range correlations. The structure factor of 
density fluctuations shows a maximum at finite wavenumber, shifting in time 
to smaller wavenumbers and corresponding to a growing correlation length. 
It agrees well with two-dimensional molecular dynamics simulations. 
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Typeset using REVTgX 



In most studies of rapid granular flows, also called the granular gas regime [J]], the inelas- 
ticity of granular collisions is assumed to be the most important feature that distinguishes 
these flows from usual liquid or gas flows. The dynamics is modeled by a single inelasticity 
parameter e = 1 — a 2 , where a is the coefficient of normal restitution. As a consequence a 
granular flow can only be maintained in driven systems, where energy is put into the system 
e.g. by gravity, shear or in vibrated layers Also quite some work has been done on the 
freely evolving granular fluid which has been shown to be linearly unstable (onset of 

clustering instability) with respect to spatial fluctuations in density, 5n(r,t) = n(r,t) — (n) 
0. In Ref. J7|, an analytic description has been given of the buildup of equal time spatial 
correlations in the flow field, 

G a p(r,t) = ^J dr'(u a (r + r',t)up(r',t)), (1) 

of a system initialized in a spatially homogeneous state. The theory is based on fluctuating 
hydrodynamics and the assumption of incompressible flow, V • u = 0. This theory yields 
predictions, including long range tails ~ r~ d , that for nearly elastic particles (e < 0.2) agree 
well with 2-D molecular dynamic simulations up to large distances. Here we will extend 
the theory to the general, compressible case, allowing us to calculate longitudinal velocity 
and density fluctuations. In fact, the structure factor corresponding to density fluctuations, 
S nn (k,t) = V~ l (5n(k, t)Sn(— k, t)), has been analyzed before by Deltour and Barrat [|J. 
The difference between compressible and incompressible flow is best appreciated in Fourier 
space, where velocity correlations are described by the tensor S a p(\i, t). Both G a p(r, t) and 
its Fourier transform S a p(k, t) = V -1 (w a (k, t)up(— k, t)) are isotropic tensors and can be 
decomposed into two scalar isotropic functions in the following way: 

G a/3 (r, t) = r a rpG\\ (r, t) + (S af3 - r a rp)G x (r, t) 

S a/3 (k,t) = k a kpS\\(k,t) + (5 a p - k a kp)S±(k,t), (2) 

where carets denote unit vectors. In a system of elastic hard spheres (EHS) for times 
larger than the mean free time to, the correlation functions are given by the equilibrium 
values, i.e. G a p(r,t) = [T/mn]6 a p5(r), containing self-correlations only, and G nn (r,t) = 
n5(r) +n 2 [g(r) — 1], where g(r) is the pair distribution function in thermal equilibrium. For 
convenience we substract self-correlations and introduce the functions G^Jr, t) = G a p(r, t) — 
[T(t)/mn\8 a p8(r) and S^(k,t) = S a p(k,t) — [T(t)/mn]S a p. Note that T(t) is measured in 
energy units (ks = 1). The structure factor of transverse velocity fluctuations, S±(k,t), was 
calculated analytically in Ref. M and shown to yield a long range r~ d tail in G±(r,t) and 
G\\(r,t) in case the fluctuations in the flow field are incompressible, i.e. St(k,t) = 0. 

In this Letter the structure factors S a p(k,t) and S nn (k,t) and corresponding spatial cor- 
relation functions G Q/ g(r, t) and G nn (r, t) will be calculated and compared with 2-D molecular 
dynamics simulations for inelastic hard disk systems. We show in particular by explicit calcu- 
lation that for small inelasticity (e ^ 0.2) S^(k, t) is essentially vanishing for all wavenumbers 
except at very small k values (k ^ 1/£||) ; where the assumption of incompressible u fluctu- 
ations, made in Ref. [0], breaks down. Consequently, the most important role of S^(k,t) is 
to provide an exponential cutoff for the r~ d tail at the largest scales r > 27r£ii. At larger 
inelasticities the contributions from St(k,t) modify G\\(r,t) and G±(r, t) significantly at all 
distances. 
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The hydrodynamic equations for the unforced inelastic hard sphere (IHS) fluid possess 
an exact solution, the homogeneous cooling state (HCS), with a homogeneous density n and 
temperature T(t), and vanishing flow field. Energy is dissipated at a rate 2 r yQ(jT, where 
7o = e/2d and where the collision frequency u = Q d x(n)no- d ~ 1 Jt~/ Tim is calculated from 
the Enskog-Boltzmann equation [|J for a dense system of hard disks or spheres (d = 2,3). 
Here Q d = 2n d / 2 /T(d/2) is the surface area of a d- dimensional unit sphere, a and m the 
sphere diameter and mass, both of which are set equal to unity, and x{n) the pair correlation 
function at contact. For detailed definitions and derivations we refer to ||. In the following, 
we will assume that IHS hydrodynamics can be described by the standard hydrodynamic 
equations supplemented by an energy sink term, which is evaluated in the local homogeneous 
cooling state. The equations of change for the macroscopic fields then become 

d t n + V • (nu) = 

9 f u + u-Vu = --V-n (3) 

n 

d t T + u-VT = -- ^-(V J + n : Vu) - 2 7o cu[n, T]T, 
dn 

where all fields and fluxes depend on (r, t). A possible justification of these equations to 
lowest order in e can be found in Ref. ||, as well as a discussion of higher order terms. The 
pressure tensor II = IIo + IIi is given by n = pi = nT(l + n d xna d /2d% with I the identity 
matrix, and IIi = — 2i]{ Vu} s — £(V • u)I, with {Vu}^^ = \{V a up + Vpu a — 2(5 a/ gV -u/d), 
and the heat flow J = — kVT. Here 77, ( and k are, respectively, the shear viscosity, bulk 
viscosity and heat conductivity, given by the Enskog theory for EHS with temperature 
T(t) still depending explicitly on time to account for the homogeneous cooling M. The 
equations of change for the mesoscopic fields are obtained from the above equations by adding 
fluctuating terms to the pressure tensor and heat flow, denoted by II and J respectively 
|T0| . They are characterized by a vanishing average and correlations which are local in space 
and time, the strength of which is assumed to be determined by the standard fluctuation- 
dissipation theorem: 

(n^r^fL^tr',^)) = 2T[r](6 aj 6p S + S aS 8 Pj ) 

+(C - ^)<WM<K r " r ')^ - 
(j a (r,t)jp(r',t')) = 2KT 2 5 a p5(r-r')5(t-t'), (4) 

with transport coefficients depending on T(t). We are interested in the buildup of cor- 
relations between spatial fluctuations in a system that is prepared in a homogeneous 
state at an initial temperature T and reaches the HCS within a few mean free times 
to = l/w[n,To\. Therefore, we can linearize the above equations around a homogeneous 
density n and a temperature T(t) = T /[l + r yot/t } 2 , and a vanishing flow field. At this 
point it is convenient to make the change of variables, dr = ou[n,T(t)]dt, where r is the 
average number of collisions a particle has suffered within a time t, 5n(r,t) = ndu(r,r), 

u(r,t) = y / T^)w(r,r), 5T(r,t) = T(t)69(r,T), tl{r,t) = nw[n,T(t)]^T(t)7r(r,r) and 

J(r, t) = nuj[n,T(t)]T(t)j(r,T). In these new variables the noise strengths of the reduced 
fluctuating pressure tensor n and heat flow j are time independent, and the equations 
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of change for the mesoscopic Fourier modes 5u(k, r), w(k, r) and 58(k, r) become or- 
dinary differential equations with time independent coefficients (valid for Mq ^5 1 where 



l = j2T(t)/oj[n,T(t)\ is the time independent mean free path): 

dSu iklo 



dr ~ ^2 

dw\ n 



U'l 



7o(l - k £ ± )w ±a - ikn a i 
5v — ikifu 



iklo 



\/2 \nTxT / 



856 



- M1 + k^e- l ^(^L) m 



dr ' uv M/ VrfnTy 

-2 lo (l + -%5u-tkh. (5) 

Here we have introduced the time independent correlation lengths £z and £t, defined 
by £j_ = v/oj^q with v = r]/mn the kinematic viscosity, £f = [2z/(c/ — l)/c? + (/mn}/uj'~fQ 
and £f> = 2K/dnwy , and the isothermal compressibility = (dn/dp) T /n. The subscript 
a in the equation for w_|_ refers to any of the {d — 1) directions perpendicular to k, and the 
subscript I denotes the longitudinal direction along k. To calculate the structure factors we 
also need the Fourier modes with k replaced by — k. 

Since the transverse velocity is decoupled from the other modes, its structure factor 
S±(k,t) = (u± a (k, t)u± a (— k, t))/V can be obtained in the analytic form |7j] 

which is valid for kl < 1. The same result has been obtained from a more microscopic 



approach, using ring kinetic theory [11 



The density, longitudinal velocity and temperature modes are coupled and their equations 
of change can be written in matrix representation as 

^(k)=M(k)^(k) + f(k), (7) 

where ip is the column vector with components ipi = 5u, ip2 = Wi and ips = 56, and the 
hydrodynamic matrix M and the noise vector f are given by Eqs. (|5|). Note that the elements 
M 31 (k) and M 33 (k), entering the temperature equation, depend on the energy dissipation 
term. In this notation the equal time correlations obey the equation of change 

^-(^(k,r)^(-k,r)) = Af Q7 (k)(V 7 (k > r)V /3 (-k,r)> 
+Af^(-k)^ (k,r)V 7 (-k,r)> + C a/3 (fc), (8) 

where a, [3, . . . — 1, 2, 3 label the components 5v, wi and 56. These equations constitute a set 
of 3 x 3 linear ordinary differential equations, of which only 6 are independent. The matrix of 
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noise strengths C al 3(k), defined through ( j a (k, r) fp(— k, r')) = C a( g(&)5(T — r'), has only two 
nonvanishing components, namely C22 = 2V^7ofc 2 £ 2 /n and C33 = AV'yok 2 ^/ 'dn. We have 
solved the above set of equations numerically, starting from initial equilibrium correlations, 
of which the only nonvanishing ones are (ipi(k, 0)ipi(— k, 0)) = VTxt, (^(k, 0)tp2(— k, 0)) = 
V/n and (03 (k, O)03(— k, 0)) = 2V/dn (for A; 7^ 0). The most important new results with 
respect to Ref. J7J] are the structure factors S\\(k, t) and S nn (k, t), and the correlation function 
G nn (r, t). In Fig. 1, we show the results for these structure factors, including Sj_(k,t), for a 
system with area fraction = 0.245 and a = 0.9 together with the results from a molecular 
dynamics simulation of 50000 inelastic hard disks. We observe that S\\(k — > 0, t) = S±(k — > 
0,t), implying for large distances an asymptotic behavior G Qj a(r, t) ~ Sj_(k — > 0,t)5 a /35(r), 
and thus the absence of algebraic long range correlations on the largest scales (r ^> 27r£||). 
Therefore, we can already conclude that the asymptotic behavior of Gj_(r, t) and G\\(r,t) 
cannot be r~ d . Instead the r~ d tail obtained in Ref. [|7| describes intermediate behavior 
which is exponentially cut off at a distance determined by the width of S^(k, t). This width 
can be estimated from the eigenvalues of the hydrodynamic matrix, more precisely from the 
dispersion relation of the heat mode ||, which is a pure longitudinal velocity wi for k — > 0. 
To second order in k its dispersion relation is given by zn{k) = 7o(l — k 2 £, 2 ), with 



I 2 

k 2 7o 2 




P 



uTxt n T \ xdn dnT 



(9) 



Note that £m ~ 1/e for small inelasticity, whereas £j_ ~ £j ~ £t ~ [ see ^S- 2(c)]. 

To a good approximation S'||(/c,i) for small wavenumbers is given by expression (|5|) with 

replaced by This approximation is excellent up to wavenumbers where the exact 
numerical result for S\\(k,t) shows a little dip (see Fig. 1, r = 19.4, k ~ 0.1). At about 
the same wavenumber the structure factor S nn (k,t) reaches its maximal value, which grows 
in time. The exact position of this maximum shifts in time to smaller wavenumbers corre- 
sponding to a growing correlation length. This can be explained by the following argument: 
for k — > density fluctuations Sn(k, t) are decoupled from the heat mode and we expect that 
Snn{k — > 0,t) remains at its initial equilibrium value; at small, but nonvanishing k, Sn(k, t) 
couples in O(k) to the unstable heat mode and the maximum of fcexp[2z^(fc)r] shifts in 
time to smaller wavenumbers. 

The estimate for S nn (k, t) ~ S nn (k, 0) exp[2^(/c)r], used in Ref. 0], differs in two aspects 
from our predictions: (i) it neglects the wavenumber dependence of the coupling of density 
fluctuations to the heat mode, giving for S nn (k,t) a decreasing function of k, and therefore 
cannot explain the growing correlation length; (ii) it neglects the fluctuating parts of the 
pressure tensor and heat flow, causing only numerical deviations from our prediction. Note 
that seven out of the eight sets of data points shown in Fig. 9 of Ref. J|] are in the crossover 
or nonlinear time regime, which is estimated in Ref. [0] to occur at r cr ~ 65 for a = 0.9 and 
= 0.4 and where our linear theory breaks down. 

Using the above approximation for S\\(k,t), the structure factor S^ik, t) can be written 

as 



n Jo 

+ (<W - k a kp) exp(-s'k 2 e ± )] , (10) 



S+ a (k,t) w fds'exp^') \k a k exp(-s'k 2 ^ 
' n Jo L 



5 



where s = 270T. If the system is thermodynamically large (L ^> 27r£||), Gt(r,t) = 
fafpG^pir^t) and G±(r,t) = (5 a p — ^a^^G^^r, t)/(d — 1) can be obtained by perform- 
ing integrals over k space; the resulting Gt(r,t) and G±(r,t) can then be expressed as 
integrals over simple functions. Here we only quote the results for d = 2. Using 

,/C| sin 2 ^e iqx - S(?2 = -^[1 -exp(-x 2 /4s)], (11) 



(2tt) 2 2ttx 
where cos 9 = q • x, we obtain 



?? 




m A f s , , s > I x \\ \ I X"\ 

+ 2^l dSe 6XP - eXP "47 



(12) 



for A =||, _L, where x\ = r/£\, m\\ = 1 and m_|_ = —1. The approximation of incompressible 
fluid flow of Ref. |7j is obtained in the limit £11 — > 00. It is consistent with the thermodynamic 
concept of incompressibility, i.e. xt — in @, implying an infinite speed of sound. Fig. 
2 shows g\\(x x ,s) = [n£ 2 ± /T(t)]G^(r,t) and g ± (x ± , s) = [n£ 2 jT{t)]G\{r } t) in the above 
approximation for different ratios at s = 2. In order to see r~ d behavior the second 

term of flT2] ) should dominate over the first. For £ 2 ^> £j_ , exp(— x 2 L /4s / ) can be neglected with 
respect to exp(— x 2 /4s'), and the second term behaves algebraically if x 2 ^ 4s. Restricting 
ourselves to 1 ^ s ^ 10 (i.e. moderate times where the correlations have grown above noise 
level), we can estimate that the range of algebraic decay is restricted to r < 27r£n, where the 
factor 2tt is chosen for convenience. For r > 27r£ii, the remaining exponent cuts off the r~ d 
tail. 

The predicted spatial velocity correlations G\\(r,t) and G±(r,t) have been obtained by 
numerically performing inverse Bessel transformations on the numerical results for S\\ (k, t) 
and S±(k,t). The result for G\\(r,t) corresponding to Fig. 1 includes an intermediate r~ 2 
tail and is shown in Fig. 3(a). Fig. 3(b) shows the corresponding spatial density correlation 
G nn { r ,t) obtained numerically from S nn (k,t). It confirms that the present theory correctly 
predicts the buildup of density correlations, including a negative correlation centered around 
a distance which grows in time as y/r. 

At small inelasticity (e < 0.2) the functions G\\(r,t) and G±(r,t), calculated here from 
the full set of hydrodynamic equations, differ for r < 27r£|| only slightly from the results for 
incompressible flow fields (see discussion in Ref. [0]). However, the algebraic tails ~ r~ d in 
G\\(r, t) and G±(r, t), derived in Ref. for r > 27r£j_, are exponentially cut off for r > 27r£||. 
As the correlation lengths £j_ ~ an d £11 ~ l/ e are we ^ separated for small e [see Fig. 

2(c)], there is an intermediate range of r values where the algebraic tail ~ r~ d in G»{r,t) 
can be observed. 

At higher inelasticity £y and are not well separated and, as a consequence, there does 
not exist a spatial regime in which the longitudinal fluctuations in the flow field can be 
neglected and the regime of validity of the incompressible theory of Ref. has shrunk to 
zero. Fig. 3(c) compares results from incompressible and compressible fluctuating hydrody- 
namics with simulation data for G±(r, t) at a — 0.6 and = 0.4 and confirms the validity of 
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the compressible fluctuating hydrodynamics description up to reasonably large inelasticities. 
Note that G nn (r,t) can only be calculated from the compressible theory. 
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support of the foundation 'Fundamenteel Onderzoek der Materie (FOM)', which is financially 
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FIGURES 



FIG. 1. Theoretical predictions (solid lines) for S±(k,t), S\\(k,t) and S nn (k,t) versus ka for 
4> = 0.245 (Z ^ 0.8) and a = 0.9, where £± = 4 and £|| = 17, at r = 19.4, 40 and 48.4, 
compared with results from a single molecular dynamics run of 50000 particles, implying a minimal 
wavenumber k m \ n a = 2-ira/L ~ 0.016. All structure superimposed on the plateau values presents 
long range correlations of dynamic origin; equilibrium structure in S nn is only present for ko^lix. 

FIG. 2. (a) 10 log <7|| (x±, s = 2) versus 10 logxj_ from incompressible fluctuating hydrodynamics 
(solid line) and present theory (|i~2|) (dashed lines from left to right: £||/£± = 1,2,5,10); as £n/£± 
decreases, the r~ 2 tail is cut off exponentially at smaller distances and finally disappears at £n = 
(b) g±(xj_,s = 2) versus x± = r/^±; the depth of the minimum decreases with decreasing £y/£± 
and finally disappears at £ii = (c) versus area fraction <p. 

FIG. 3. (a) 10 log [| C|| \/T] versus 10 logr. (b) G nn versus r; the same parameters a,0,r as in 
Fig. 1 are used for (a) and (b). (c) G±/T versus r for <j) = 0.4 (l ~ 0.34), a = 0.6 (£± = 1.46, 
^11 = 3.8) at r = 20, 40 and 60; in (a) and (c) the solid (dashed) line is the prediction from 
compressible (incompressible) fluctuating hydrodynamics. 
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